KIVA Family of Codes:
A General Introduction

by 
Norman L. Johnson
Los Alamos National Laboratory
(norman@santafe.edu)

What is KIVA?
 A KIVA is a cylindrical ceremonial cave used by the Pueblo Indians in the Southwest United States.  It’s also the most recent family of computer codes for simulating reactive flows typically found in combustion engines, developed by the Fluid Dynamics Group (T-3) in the Theoretical Division at Los Alamos National Laboratory in New Mexico.  While many combustion codes treat reactive flows, KIVA is somewhat unique in that it includes the capability to model liquid sprays, typical of fuel injection systems.  This introduction to the KIVA family of codes begins with a summary of the physics and relevant concepts of computational fluid dynamics in order to establish a common vocabulary.  The lineage of KIVA is reviewed.  And finally the approach to the gas dynamics is presented.  
Problem Definition: Governing equations, constitutive equations and parameters, boundary conditions, and numerical parameters
A simulation of a combustion engine includes a wide variety of physical phenomena, all typically interconnected.  In addition to the physics of a fluid dynamics code, the added complications of species diffusion, chemical kinetics, heat transfer, multiphase flow, phase changes and others are possible.  Modeling of physical phenomena requires the determination of the governing sets of equations (such as the evolution equations for mass, momentum and energy).  To apply the governing equations, it’s also necessary to develop constitutive relationships, which typically relate the dynamic and thermodynamic state to forces or fluxes in the governing equations (such as an equation of state, Newton’s law of viscosity, Fourier’s law of heat conduction).  Constitutive equations require material parameters, obtained by theory or by an experimental measurement.  The final requirement of a problem definition is the extent of the domain of a problem and the specification of the boundary conditions.  Once the governing equations, constitutive equations with material parameters, and the appropriate boundary conditions are determined, these equations, typically describing a continuum, must be recast in a form that can be solved on a computer. 
Even if you are not concerned with the modification of the governing and constitutive equations, it is essential to know the assumptions or limits of their applicability in order to apply the model.   Excluding numerical difficulties, violating these assumptions and limits will be the most likely source of inaccuracies in a model.  
Common Concepts to a Fluid Dynamics Code
To establish a common vocabulary in order to discuss the KIVA family of codes, an introduction of the basic terms used in computational fluid dynamics (CFD) codes follows.  A more complete discussion can be found in most CFD texts, such as “Computational Techniques for Fluid Dynamics” by C. A. J. Fletcher (Springer-Verlag, 1991, 2nd Edition.). 
Difference equations from continuum equations
The governing equations of the dynamics of the an engine system are almost always differential or integral equations describing a continuum.   Some methodology is required to transform the governing equations to discrete equations, as required for a computational solution.  The resulting discrete equations describe the time and spatial evolution of the system at a collection of computational points.  The methodology typically used in the KIVA codes is the finite volume method, a more modern version of the finite difference method.  The finite volume method spatially decomposes the domain into elements (or cells, in the typical Los Alamos nomenclature), transforms the volume integrals into surface integrals and then replaces the continuous surface contributions with discrete contributions as determined from the computational points.   The resulting set of discrete equations are not unique for a given methodology; for example, with the proper choice of basis functions and order of integration, the finite element method can usually reproduce the same discrete equations as a finite volume method.  
Meaning of computational points
All CFD codes have spatial points where information is located.  One major distinction among various CFD codes is whether the computational points represent a fixed mass that moves with the material - a Lagrangian representation, or a fixed position in space with mass moving through the point - an Eulerian representation.  Numerical methods based on Lagrangian points have the advantage of no diffusion of information due to material movement, but most implementations suffer from mesh entanglement in shearing flows.  Methods based on the Eulerian representation have the advantage of no mesh distortion, but can suffer from an accurate description of interfaces and diffusion of information from mass transport.  An alternative to either of these is the Arbitrary Lagrangian-Eulerian (ALE) representation.  Here, a computational point can have an arbitrary velocity.  Depending on the choice of the velocity, the point can be Lagrangian, Eulerian or a specific choice as needed to optimize some desirable mesh configuration, such as reduced mesh distortion or a more resolved mesh, capturing a feature of interest.  
A common confusion that arises for computational methods is whether a computational point represents a value of the continuum field at the location of the point or represents some volume-averaged value.  This confusion is particularly true for the less rigorous finite difference method.  Sometimes it is necessary to have different interpretations at different times in the numerical treatment; this is particularly true for computational points along a boundary.  In KIVA, the volume-averaged values are the standard volume interpretation of the points, with an occasional point interpretation as necessary. 
Location of dynamic and state variables
Associated with each computational point is information required for the numerical method.  Where information is located can determine many aspects of the method - the ease of implementation, the accuracy, the degree of conservation.  For the KIVA family of codes, all information is staggered on the mesh, meaning that the velocities are located separately from the thermodynamics state (internal energy, pressure, density, etc.).  State information (pressure, internal energy, density) is located at cell centers; velocities are located at cell edges.  Alternatively, all variables can be located at the same computational point and then are said to be co-located.  Staggered numerical methods have the typical advantage of tightly coupled momentum and pressure fields but can suffer from lack of conservation, due to the control volumes for the momentum and energy not being spatially coincident. 
Types of meshes and connectivity
In most CFD codes the connections between the computational points define broad classes of codes and strongly determine the ease or complexity of the implementation of the numerical method or the construction of a mesh for a complex geometry.  A mesh where the connections are unchanged from the initial time is said to have fixed connectivity, in contract to variable connectivity codes in which the connections change continuously or at certain events in the simulation (e.g., opening of a port).  If a mesh is viewed with only its connectivity considered, ignoring the spatial relationship of the computational points, then it is called a logical mesh.  The figure below (missing in this document) illustrates how one logical mesh can represent widely different spatial meshes. If the connections in a logical mesh follow a simple pattern throughout the mesh, then the mesh is called structured.  The connections in an unstructured mesh have no pattern.  Typically CFD codes based on unstructured meshes have more flexibility in mesh generation, but more complexity in the code structure.  Because structured meshes have well defined paths through the logical mesh, logical coordinates can be defined along the mesh directions and typically have integer values or an integer plus some fraction.  These logical coordinates have no relationship to the spatial coordinates.  KIVA-2 uses a single logical block of mesh, comprised of hexahedrons.  KIVA-3 uses a collection of logical blocks of mesh of hexahedrons in which the connectivity can change during the simulation. 

	Logical mesh	Spatial mesh #1	Spatial mesh #2
Approaches and properties of numerical solutions
Discrete equations can be cast and solved in different ways that are compromises between stability, accuracy, consistency and efficiency of the solution.  
	•	Stability refers to the property of a solution to remain bounded.  
	•	Accuracy refers to how close an approximate numerical solution approaches an analytical solution of the discrete equations.  
	•	Consistency refers to how close the analytical solution of the discrete equations duplicates the solution to the governing equations. 
	•	Efficiency generally refers to the computational ease with which the numerical equations are solved.  
In solving for the time evolution of a problem, if the future state can be solved directly from the current state, then the method is said to be explicit.  A time step that advances the evolution of the problem is typically chosen based on the stability of the solution.  Alternatively, if the solution of the future state also requires information at the future state, then the method is said to be implicit and typically requires a matrix inversion solutions or an iterative solution.  The time steps in implicit solutions are typically larger than for explicit solutions and are based on accuracy considerations.  Implicit solutions also are more computationally intensive and are therefore less efficient for the same size of time step . 
A test for the accuracy of a solution (but not consistency!) is to reduce the time step or increase the mesh resolution.  If the two solutions are not significantly different, then the solution is considered accurate with respect to the time step or mesh size.  This test for accuracy may not work if the solution includes multiple time or length scales.  
Family Lineage of KIVA
The development and complexity of KIVA has its origins in a long lineage of computer codes developed at Los Alamos for simulating reactive flows.  Each new code or version has added additional numerical or physical capability, with the goal that all new capabilities are integrated with existing features.  The following is a list of the past and future lineage of KIVA codes. 
RICE  (1975)
A two-dimensional Eulerian reactive flow code with turbulence modeled as an eddy diffusivity. It included kinetics and was partially implicit.
APACHE (1979)
The RICE code with arbitrary shaped cells. 
CONCHAS (1979)
An ALE version or APACHE with a turbulent subgrid scale model.
CONCHAS-SPRAY (1982)
The CONCHAS code with a statistical treatment of sprays, wall treatments for turbulent heat and momentum transport, kinetic and equilibrium chemistry.
KIVA (1985)
The CONCHAS-SPRAY with extensions to fully three-dimensional flows and droplet collision and breakup.  The option of the elimination of acoustic waves by subcycling was added. 
KIVA-II (1989)
The second version of KIVA significantly extended the numerics and submodels of KIVA by the addition of an implicit treatment of diffusion and pressure, k-e model for turbulence with advection of length scale instead of turbulence dissipation rate, a material library for gases and fuels, more boundary conditions, and a model for aerodynamic breakup of sprays. Although the mesh was limited to a single hexahedron of logical mesh, the capability of deactivating cells allowed for complex geometries.  
KIVA-3  (1993, publicly released version)
KIVA-3 contains the same numerics and physics of KIVA-II, but the single logical block of mesh of KIVA-II is replaced by a mesh of interconnected logical blocks of mesh with variable connectivity, thereby enabling an improved piston treatment, the addition of ports and more efficient use of mesh than in KIVA-2.  Both a workstation and Cray version of KIVA-3 were developed.  (Note that the Cray version of KIVA-3 is currently only available to government facilities, due to a cooperative research agreement with Cray Research.  The workstation version is still available.) Companion codes of a mesh generator and graphics post-processing code accompanied the release of KIVA-3.  
KIVA-3  (next released version - Summer 1996)
A new version of KIVA-3 will be released when the supporting documentation has been written, estimated to be in early summer 1996.  To eliminate the maintenance of two versions of the code, in the future only the workstation version will be released; the Cray version will no longer be supported.  The workstation version works on all Crays, with a small penalty in run time. The major addition in capability of the new KIVA-3 version is a valve model based on Hessel’s implementation (UW Madison Ph.D. thesis, 1993).   Because the approach resolves the valves without additional approximations, the accuracy of the valve treatment is of comparable accuracy to the flow dynamics in the rest of the code.  The valve capability allows canted valves (the valves need not be aligned with the piston axis), with individual orientations and lift histories and automatically handles spatially overlapping valves and pistons.  Mass fluxes through the valves are calculated and output during the run.  Also included with the KIVA-3 release will be the accompanying pre- and post-processors that support the valve capability.  Because of the implementation of the valve model, the interior mesh in the neighborhood of the valve has certain restrictions, which are supported in the mesh generation in K3PREP. All previous KIVA-3 capabilities have been retained and are compatible with the valve model. 
Other new features of the release are the inclusion of a mixing-controlled combustion model and a soot production model.  A new sorting subroutine is included that is linear in the number of nodes and preserves the original storage sequence. The law-of-the-wall heat transfer model has been modified and is in closer agreement with experiments.  A better description of the temperature dependence of surface tension and new fuel options have been added.  An option has been included that improves the performance in low Mach number flows.  And finally a no-hydro option has been added that runs through the simulation without fluid flow, which allows a rapid determination of mesh integrity during port, piston and valve movement. 
KIVA-F90 (Public release awaiting documentation)
Parallel with the effort to continue the maturation of KIVA-3, future versions of KIVA are being developed at Los Alamos.  These codes use largely the same numerics as KIVA-II but address the requirements of parallel computer architectures and requirements of modern mesh generation codes.  KIVA-F90 is a complete rewrite of KIVA-II using FORTRAN 90 and executes on workstations, massively parallel architectures, and supercomputers without modification.   The use of FORTRAN 90 significantly simplifies the coding and results in an easier to maintain code.  The code is also highly modular, with equation of state and boundary conditions in separate subroutines.  Message passing is not part of the code but is easily implemented due to the structure of the code. 
KIVA-AC (Under development)
KIVA-AC, now under development, is an unstructured mesh version of KIVA-F90 that will support combinations of tetrahedrons and hexahedrons. 
Other versions of KIVA
The Engine Research Center (ERC) at the University of Wisconsin, Madison has significantly extended the physics models in KIVA-II; these are summarized in the next section.  The ERC has also developed companion codes to aid in problem setup, mesh generation, and post-processing. Information on the availability of these should be directed towards Prof. Reitz (reitz@engr.wisc.edu).  
A parallel, optimized KIVA-3 for Silicon Graphics platforms is also available from Energy Science and Technology Software Center,  (615) 576-2606, Fax:  (615) 576-2865 or estsc@adonis.osti.gov.  The coding has been optimized for R8000 based machines: POWER CHALLENGE and POWER Indigo 2, but will execute on any SGI platform.  The code will execute on parallel processors if available.  
Comments on the stochastic nature of KIVA
Typical regimes of the solution of the governing equations and numerical implementation in KIVA result in deterministic solutions, meaning that with the identical initial problem definition, the same solution will evolve.  There are two exceptions to the deterministic solutions.  One is when the flow is chaotic, and small changes in the initial conditions or precision of the solution (such as a change of computing platform) can lead to appreciable changes in the final results.  These types of flows have been well studied and generally do not arise in typical KIVA engine simulations due to the diffusive nature of the turbulent model.  The other non-deterministic aspect of KIVA is the stochastic treatment of the spray.  The solution to the stochastic equations for the drop dynamics and the nature of the implementation of the equations result in an increased sensitivity to the slightest variation in the initial conditions of the setup, the precision of the computer or the choice of a seed for the random numbers.  As the number of parcels is increased in the representation of the spray, this sensitivity is reduced but convergence to a machine-independent or seed-independent solution can be quite slow - some fractional power of the number of parcels.  This stochastic nature of KIVA is not generally appreciated, and the sensitivity and convergence of the simulation results should be determined by increasing the number of parcels in the simulation, by small changes in the spray setup (such as the spray location) and different seeds for the random number generator.  
KIVA Solution Approach
Typical of most Eulerian and ALE codes, including KIVA, is the division of the solution of the governing equations into a Lagrangian phase and a remap or rezone phase.  This greatly simplifies the discrete equations that must be solved.  The Lagrangian phase fixes the mass of each element and solves for the new state, allowing the mesh to move with the gas.  The implicit approach in KIVA in the Lagrangian phase requires multiple inner iteration loops on specific physics and a larger iteration loop to coordinate all of the physical effects.  If the time step is too large or the test for convergence inappropriate, any of these iteration loops can fail to converge.  In practice, the temperature iteration is the most likely to exhibit problems.  (For possible solutions, see the next section.)     The remap phase holds the mass fixed, moves the mesh to the desired position and exchanges (or fluxes) the mass between the elements in order to account for the mesh movement.  Subcycling (the breaking up of the time step into a number of smaller time steps) of the remap phase is usually required because of the typically large time steps used in the implicit Lagrangian phase.  
KIVA Gas Dynamics
The implementation of the gas dynamics is fairly standard for a reactive flow simulation code, with the exception of the additional considerations for the interaction of the spray model with bulk flow and turbulence models.  A detailed discussion of the gas dynamics and the numerical implementation can be found in the KIVA-II manual (KIVA-II: A computer program for chemically reactive flows with sprays by A. A. Amsden, P. J. O’Rourke, and T. D. Butler. Report LA-11560-MS, Los Alamos National Laboratory, 1989).  After a brief description of the implementation of the gas dynamics in KIVA, the remainder of the discussion focuses on the parameters controlling the gas dynamics and general recommendations.  
The gas dynamics are described with the continuum equations for mass, momentum and energy, with options for either a laminar viscosity, a single turbulent transport equation (the subgrid scale (SGS) model) or the modified k-ε transport equations.  In the presence of a spray, the mass, momentum and energy equations of the gas contain additional source terms that account for the presence of the spray.  The major assumptions of the formulation are (1) state relations are those of an ideal gas mixture, (2) the laminar viscosity varies only with temperature, (3) Fick’s law of diffusion applies with a single diffusion coefficient for all species and (4) Fourier’s law of heat conduction applies with a constant value of the thermal conductivity.  
The following is a list of the input parameters affecting the gas dynamics with comments on their usage.  Consult the KIVA-II  and KIVA-3 manual ( KIVA-3: A KIVA Program with Block-Structured Mesh for Complex Geometries by A. A.  Amsden.  Report LA-12503-MS, Los Alamos National Laboratory, 1993) for a more complete description of the parameters.
PGSSW - This parameter activates the Pressure Gradient Scaling (PGS) method for better efficiency in computing low speed flows in which acoustic waves are not important and the system does not have applied pressure boundaries.  The code will automatically disable the PGS method if an applied pressure boundary is used.  Use of the PGS method increases the size of the possible time step and improves the convergence efficiency during the pressure iteration. 
EPSY, EPSV, EPSP, EPST, EPSK, EPSE - These parameters are the allowed relative errors in the implicit solution of the mass, momentum, pressure, temperature, turbulent kinetic energy and turbulent dissipation.  If these are too small, then the code will not converge in the maximum allowable iterations.  Typically these are not modified.  If simulation repeatedly exceeds the maximum number of iterations, then these can be modified or the form of the test in which they are used should be modified, but only after other options are considered.  Often a failure of converge is cause by a rapid change in the flow dynamics, and the past time step is too large for the new situation.  First try to reduce the time step and see if the problem will now converge.  The addition of node coupling can also help (see below).  In some problems, better convergence can also be obtained by making the solution more implicit (see the KIVA-3 manual for details). 
PARDON, A0, B0 - These parameters select and modify the different fluxing options in the rezone phase.  Generally these are left unmodified and are defaulted to values that emphasize stability, rather than accuracy.  Modification of these parameters should be considered if the simulation involves the transport of sharp gradients through the mesh. 
ANC4 - This parameter determines the degree of node coupling as required to dampen pressure and velocity oscillations in a simulation.  Due to the improved stability of the current numerics in KIVA, the node coupling should not be used because of its strong diffusive nature.  The addition of node coupling can also help simulations that are exceeding the maximum iterations before converging. 
ANGMOM - A classic shortcoming of Eulerian codes is their inability to sustain a swirling flow, in comparison to Lagrangian codes, due to the lack of conservation of angular momentum during advection.  The KIVA code has the option of conserving angular momentum about a central z axis at the expense of linear momentum conservation.  This option should be used with care: it is accurate only for the idealized swirling flow that is centered on around the z-axis.  The option may unrealistically alter any other type of swirling or tumble flow.  It also is most accurate for meshes with an axis along the z direction (a polar or “pie-shaped” mesh as in the earlier figure labeled “Spatial Mesh #2”).  
TKESW, SGSL, TKEI, UNISCAL, AIRMU1, AIRMU2 -  These parameters select between the different viscous and turbulent options and initialize the turbulent fields.  The use of the turbulent models is beyond the scope of this text, but some general comments can be made.   The intent of the turbulent models in KIVA is to capture the chaotic turbulent eddies in continuum equations more efficiently than resolving them with the Navier-Stokes equations.  This is accomplished at the expense of many assumptions about the nature of the turbulence (a fixed spectrum of eddy scales, isotropy, homogeneity, and others).  More often than not, these assumptions are not appropriate, but other practical alternatives are not available.  Special care must be taken when the simulation results are directly dependent on the turbulence levels, as in the case of the flame speed in turbulent combustion.  One technique for determining the reasonableness of the turbulence fields is to examine plots of the length scale, a standard K3POST plot, and to look for unrealistic values.